TD3 - 17/05/21 - Simulation et Monte Carlo

  • Exercice 7 : MCMC
  • Exercice 9 : Cross-Entropy Method

Exercice 7 : MCMC

Rappel :

image.png

7.1.

  • Quelle loi reconnaissez-vous ?

image.png

7.2.

  • Calculer les lois conditionnelles de chaque composante.

image.png

image.png

  • Mettre en oeuvre le Gibbs sampler correspondant

image.png

image.png

Code

In [12]:
from scipy import stats
import numpy as np

def rand_gauss_gibbs(n, rho):
    
    # initialisation : nous pouvons choisir autre valeur
    X = np.zeros((n, 2))
    
    # itération gibbs sampler
    for t in range(n-1):
        X[t, 0] = stats.norm.rvs(rho*X[t, 1], np.sqrt(1-rho**2), 1)
        X[t+1, 1] = stats.norm.rvs(rho*X[t, 0], np.sqrt(1-rho**2), 1)
    
    # ne pas oublier la dernière simulation de X1
    X[n-1, 0] = stats.norm.rvs(rho*X[n-1, 1], np.sqrt(1-rho**2), 1)
    
    return(X)

rn = rand_gauss_gibbs(10, 0.6)
In [13]:
rn
Out[13]:
array([[ 0.60746691,  0.        ],
       [-0.15985096, -0.29046284],
       [-1.00371489, -1.22978296],
       [-2.06407901, -0.56233516],
       [-0.30027895, -0.96924788],
       [ 0.28987404,  0.68579595],
       [-0.90168825, -0.80166987],
       [-0.44625301,  0.92051056],
       [ 0.80502245,  1.13155054],
       [-1.53346552,  0.39101761]])
  • Représenter dans le plan : des contours de la densité cible et l'évolution de la chaîne de Markov simulée
In [38]:
import matplotlib.pyplot as plt

# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_gibbs_contour(n, rho):
    
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    
    # contours de la densité cible
    rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
    
    x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
    pos = np.dstack((x, y))
    
    plt.contour(x, y, rv.pdf(pos))
    plt.title(f"rho = {rho}")
    plt.xlabel(f"X1")
    plt.ylabel(f"X2")

    # représentation de la chaine de Markov simulée
    plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
    return None
In [111]:
plot_gibbs_contour(n = 500, rho = 0.5)
  • Que se passe-t-il quand $\rho \to 1$ ?
In [109]:
n = 100
rhos = [0, 0.2, 0.5, 0.7, 0.9, 0.99]

plt.figure(figsize = (16, 10))

for i, rho in enumerate(rhos):
    plt.subplot(2, 3, i+1)
    plot_gibbs_contour(n = n, rho = rho)

lorsque $\rho \to 1$, $X_2 = X_1$ p.s.

rq : on peut montrer (et voir) facilement que lorsque $\rho \to -1$, $X_2 = -X_1$ p.s.

Nous utilisons les package python statsmodels pour évaluer la performance MCMC

In [96]:
from statsmodels.graphics.tsaplots import plot_acf
import pandas as pd

def plot_gibbs_ACF(n, rho, dim = 0):
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    # plot
    plot_acf(rn[:,dim])
    plt.title(f"Autocorrelation pour X{dim+1}")
    return None

def plot_gibbs_trace(n, rho, dim = 0):
    # simuler X1 et X2
    rn = rand_gauss_gibbs(n, rho)
    # plot
    plt.plot(np.arange(n), rn[:,dim])
    plt.title(f"Trace X{dim+1}")
    return None
In [179]:
n = 1000
rho = 0.5
plot_gibbs_ACF(n, rho, dim = 0)
plot_gibbs_ACF(n, rho, dim = 1)
plt.show()
In [178]:
n = 1000
rho = 0.5

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()

Que se passe-t-il si $\rho$ proche de 1 ?

In [177]:
plot_gibbs_ACF(1000, 0.99, dim = 0)
plot_gibbs_ACF(1000, 0.99, dim = 1)
plt.show()
In [176]:
n = 1000
rho = 0.99

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_gibbs_trace(n, rho, dim = 0)
plt.subplot(1,2,2)
plot_gibbs_trace(n, rho, dim = 1)
plt.show()

7.3.

  • Programmer un algorithme RWHM (Hastings-Metropolis avec marche aléatoire)
In [131]:
def RWMH(n, rho, sigma):
    X = np.zeros((n,2))

    for t in range(n-1):
        
        # proposition
        x = stats.multivariate_normal(X[t,:], sigma).rvs(1)
        
        rvn = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
        
        if np.log(np.random.rand()) < np.log(rvn.pdf(x)) - np.log(rvn.pdf(X[t,:])):
            X[t+1,:] = x
        
        else:
            X[t+1,:] = X[t,:]

    return X

n = 10
rho = 0.8
tau = 0.8
X = RWMH(n, rho, tau * np.eye(2))
X
Out[131]:
array([[ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [ 0.        ,  0.        ],
       [-1.31137076, -0.57368188],
       [-1.44696416, -0.39851567],
       [-1.44696416, -0.39851567],
       [-1.44696416, -0.39851567],
       [-1.44696416, -0.39851567]])
In [142]:
# fonction qui affiche les contours de la densité cible et l'échantillon simulé
def plot_RWMH_contour(n, rho, tau, sigma = np.eye(2)):
    
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    
    # contours de la densité cible
    rv = stats.multivariate_normal(np.zeros(2), np.array([1, rho, rho, 1]).reshape(2, 2))
    
    x, y = np.mgrid[min(rn[:, 0]):max(rn[:, 0]):.1, min(rn[:, 1]):max(rn[:, 1]):.1]
    pos = np.dstack((x, y))
    
    plt.contour(x, y, rv.pdf(pos))
    plt.title(f"rho = {rho}, tau = {tau}")
    plt.xlabel(f"X1")
    plt.ylabel(f"X2")

    # représentation de la chaine de Markov simulée
    plt.scatter(rn[:,0], rn[:,1], color = "steelblue")
    return None

def plot_RWMH_ACF(n, rho, tau, sigma = np.eye(2), dim = 0):
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    # plot
    plot_acf(rn[:,dim])
    plt.title(f"Autocorrelation pour X{dim+1} avec rho = {rho}, tau = {tau}")
    return None

def plot_RWMH_trace(n, rho, tau, sigma = np.eye(2), dim = 0):
    # simuler X1 et X2
    rn = RWMH(n, rho, tau*sigma)
    # plot
    plt.plot(np.arange(n), rn[:,dim])
    plt.title(f"Trace X{dim+1} avec rho = {rho}, tau = {tau}")
    return None
  • voici ci-dessous les paramètres que nous allons analyser :
In [271]:
taus = [0.001, 1, 5, 100]
rhos = [0.2, 0.9, 0.99]
In [272]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = 1000, rho = rhos[0], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
In [273]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = 1000, rho = rhos[1], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
In [274]:
plt.figure(figsize=(16, 4))

for i, tau in enumerate(taus):
    plt.subplot(1, 4, i+1)
    plot_RWMH_contour(n = 1000, rho = rhos[2], tau = tau)
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
In [275]:
n = 1000
plot_RWMH_ACF(n, rhos[2], taus[2], sigma = np.eye(2), dim = 0)
plot_RWMH_ACF(n, rhos[2], taus[2], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
In [276]:
n = 1000

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[2], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[2], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
In [277]:
n = 5000

plt.figure(figsize=(16, 4))

plt.subplot(1,2,1)
plot_RWMH_trace(n, rhos[2], taus[2], sigma = np.eye(2), dim = 0)
plt.subplot(1,2,2)
plot_RWMH_trace(n, rhos[2], taus[2], sigma = np.eye(2), dim = 1)
plt.show()
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':
/Users/Faugon/opt/anaconda3/envs/active/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in log
  if sys.path[0] == '':

Bilan des courses :

  • Gibbs sampling est un bon outil pour construire des chaînes de Markov

  • Cependant lorsque les v.a marginales sont très corrélées ($\rho \to 1$) nous avons des anomalies : loi simulée différente de la loi cible

  • nous pouvons améliorer la performance de la simulation en utilisant un algorithme RWHM (où $\Sigma$ de la loi auxilière $\mathcal{N}(\mu, \Sigma)$ doit être bien choisie, pas "trop grand", pas "trop petit").

Exercice 9 : Cross-Entropy

image.png

image.png

image.png

image.png

image.png

Code :

In [264]:
from numpy.random import multivariate_normal

def rosenbrock(x):
    y = 0
    for i in range(len(x)-1):
        y = y + 100 * (x[i+1]-x[i])**2 + (x[i]-1)**2
    return(y)

def optCE(fun = rosenbrock, n = 1000, d = 3, eps = 1e-8, max_iter = 1000, tau = 0.01):
    """
    this function compute the arg max of the function 'fun' by cross-entropy method
    """
    mu = np.random.rand(d)
    sigma = 10*np.eye(d)
    
    t = 0
    while True:
        t += 1
        
        # step 1 : sample Yi
        Y = multivariate_normal(mean = mu, cov = sigma, size=n)

        # step 2 : pick the top Yi that minimize the function
        samples = - np.array(list(map(fun, Y)))
        n_top = round(tau * n)
        top_ind = np.argsort(samples)[::-1][:n_top]
        Y_top = Y[top_ind, :]

        # step 3 : estimate by MLE mu and sigma
        mu = np.mean(Y_top, 0)
        sigma = np.cov(Y_top.T)
            
        # step 4 : stopping criter
        if np.max(sigma) < eps or t >= max_iter:
            print(t)
            print(mu)
            break
    return mu
In [267]:
optCE(n = 1000, d = 3)
7
[0.99999981 0.99999906 0.99999899]
Out[267]:
array([0.99999981, 0.99999906, 0.99999899])
In [268]:
optCE(n = 1000, d = 5)
9
[0.99997798 0.9999858  0.99998097 0.99998683 0.99996077]
Out[268]:
array([0.99997798, 0.9999858 , 0.99998097, 0.99998683, 0.99996077])
In [269]:
optCE(n = 1000, d = 10)
18
[-0.20302528  0.21181113 -0.02025519 -0.04240751  0.17407112 -0.06280436
  0.35965187  0.08616139  0.03233633  0.33010423]
Out[269]:
array([-0.20302528,  0.21181113, -0.02025519, -0.04240751,  0.17407112,
       -0.06280436,  0.35965187,  0.08616139,  0.03233633,  0.33010423])
In [270]:
optCE(n = 10000, d = 10)
15
[0.99999706 0.999997   0.99999753 0.99999893 1.00000209 1.00000668
 1.00000792 1.00001085 1.00000979 1.00001088]
Out[270]:
array([0.99999706, 0.999997  , 0.99999753, 0.99999893, 1.00000209,
       1.00000668, 1.00000792, 1.00001085, 1.00000979, 1.00001088])

Bilan des courses :

  • la méthode CE peut-être utilisé soit pour une estimation soit pour une optimisation
  • pour une optimisation, on doit tenir en compte plusieurs paramètres :

    • si $n$ (taille de l'échantillon) devient grand alors l'optimisation est plus précise
    • si $d$ (dimension) devient grand alors l'optimisation est plus compliqué
  • l'initialisation de $\theta_0$ est bien sûr très important, ici $\theta_0 = (\mu_0, \sigma^2_0)$ et :

    • idéalement $\mu_0$ doit être proche de la vraie valeur (sinon on choisit une valeur déterministe ou aléatoire)
    • $\sigma_0$ doit être assez grand afin de considérer toutes les possibilités